2019 Data pulled in 2022

file <- paste0(ref_path, "/PFIRS/PFIRS 2019-2020 CT pulled 2022.xlsx")
xy_old <- read_excel(file)
xy_old %>% head()
## # A tibble: 6 × 18
##   `Burn Date`         `Burn Unit`     `Acres Burned` `Total Tons`         
##   <dttm>              <chr>                    <dbl> <chr>                
## 1 2019-03-24 00:00:00 Office                    0.01 5.4300000000000001E-2
## 2 2019-12-16 00:00:00 Mammoth Bar               0.01 0.1135               
## 3 2019-04-16 00:00:00 Foresthill Loop           0.02 9.6000000000000002E-2
## 4 2019-02-08 00:00:00 Mammoth Bar               0.05 0.5675               
## 5 2019-03-13 00:00:00 Foresthill Loop           0.05 0.24                 
## 6 2019-03-27 00:00:00 Foresthill Loop           0.05 0.24                 
## # ℹ 14 more variables: `Estimated Emissions` <chr>, Agency <chr>,
## #   `Agency Unit` <chr>, Landowner <chr>, `Fuel Type` <chr>, `Burn Type` <chr>,
## #   Latitude <chr>, Longitude <chr>, `Legal Location` <chr>,
## #   `Min Elevation` <dbl>, `Max Elevation` <dbl>, County <chr>,
## #   `Air District` <chr>, `Air Basin` <chr>

Fix column names

xy_old <- xy_old %>% 
  rename(Burn_Date = `Burn Date`) %>% 
  rename(Burn_Unit = `Burn Unit`) %>% 
  rename(Acres_Burned = `Acres Burned`) %>% 
  rename(Fuel_Type = `Fuel Type`) %>% 
  rename(Total_Tons = `Total Tons`) %>% 
  rename(Burn_Type = `Burn Type`) 

Make lat and long numeric and remove unk rows

xy_old <- xy_old %>% 
  filter(!Longitude == "UNK") %>% 
  mutate(Longitude = as.numeric(Longitude)) %>% 
  mutate(Latitude = as.numeric(Latitude))

Fix dates that should be negative

xy_old <- xy_old %>% 
  mutate(Longitude = ifelse(Longitude > 0, Longitude*(-1), Longitude))

Add year column

xy_old <- xy_old %>% 
  mutate(Year = year(Burn_Date))

Filter to 2019

xy_old_2019 <- xy_old %>% 
  filter(Year == 2019)
xy_old_2019 %>% head()
## # A tibble: 6 × 19
##   Burn_Date           Burn_Unit    Acres_Burned Total_Tons `Estimated Emissions`
##   <dttm>              <chr>               <dbl> <chr>      <chr>                
## 1 2019-03-24 00:00:00 Office               0.01 5.4300000… 0.25                 
## 2 2019-12-16 00:00:00 Mammoth Bar          0.01 0.1135     0.32                 
## 3 2019-04-16 00:00:00 Foresthill …         0.02 9.6000000… 0.6                  
## 4 2019-02-08 00:00:00 Mammoth Bar          0.05 0.5675     0.32                 
## 5 2019-03-13 00:00:00 Foresthill …         0.05 0.24       0.6                  
## 6 2019-03-27 00:00:00 Foresthill …         0.05 0.24       0.6                  
## # ℹ 14 more variables: Agency <chr>, `Agency Unit` <chr>, Landowner <chr>,
## #   Fuel_Type <chr>, Burn_Type <chr>, Latitude <dbl>, Longitude <dbl>,
## #   `Legal Location` <chr>, `Min Elevation` <dbl>, `Max Elevation` <dbl>,
## #   County <chr>, `Air District` <chr>, `Air Basin` <chr>, Year <dbl>

Delete duplicates

From Jason Branz: “Typically if there are multiple records with the same date, burn name, and acres, it’s a duplicate.”

Look at duplicates

xy_old_2019 %>% 
  select(Burn_Date, Burn_Unit, Acres_Burned) %>% 
  filter(duplicated(.))
## # A tibble: 112 × 3
##    Burn_Date           Burn_Unit                          Acres_Burned
##    <dttm>              <chr>                                     <dbl>
##  1 2019-12-05 00:00:00 SPRA                                       0.25
##  2 2019-12-10 00:00:00 Unit 2                                     0.25
##  3 2019-12-10 00:00:00 SPRA                                       0.25
##  4 2019-12-11 00:00:00 Unit 4                                     0.25
##  5 2019-12-16 00:00:00 Big Chico Creek Ecological Reserve         0.25
##  6 2019-12-16 00:00:00 Big Chico Creek Ecological Reserve         0.25
##  7 2019-12-17 00:00:00 BCCER Prop Management                      0.25
##  8 2019-12-17 00:00:00 BCCER Prop Management                      0.25
##  9 2019-01-29 00:00:00 SPRA                                       1   
## 10 2019-01-29 00:00:00 SPRA                                       1   
## # ℹ 102 more rows

Remove duplicates

xy_old_2019 <- xy_old_2019 %>% 
  distinct(Burn_Date, Burn_Unit, Acres_Burned, .keep_all = T)

Make spatial

sf_old_2019 <- st_as_sf(xy_old_2019, coords = c("Longitude", "Latitude"), crs = 4326)

Export to shapefile

st_write(sf_old_2019, dsn = "shapefiles", layer = "PFIRS_2019_pulled2022", driver = "ESRI Shapefile", delete_layer = T)
## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver
## Deleting layer `PFIRS_2019_pulled2022' using driver `ESRI Shapefile'
## Writing layer `PFIRS_2019_pulled2022' to data source 
##   `shapefiles' using driver `ESRI Shapefile'
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 6: Field Burn_Dt create as date
## field, though DateTime requested.
## Writing 1100 features with 17 fields and geometry type Point.
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Message 1: One or several characters
## couldn't be converted correctly from UTF-8 to ISO-8859-1.  This warning will
## not be emitted anymore.

2019 Data pulled in 2023

xy_new <- read_excel(paste0(ref_path, "/PFIRS/PFIRS 2017-2022 pulled 2023.xlsx"))
## Warning: Expecting date in A8155 / R8155C1: got 'NA'
xy_new
## # A tibble: 10,650 × 9
##    `Burn Date`         `Burn Unit`      Agency `Acres Burned` Latitude Longitude
##    <dttm>              <chr>            <chr>           <dbl>    <dbl>     <dbl>
##  1 2017-01-02 00:00:00 PP Ridgetop Unit Cal F…           10       38.8     -123.
##  2 2017-01-03 00:00:00 Fairway/Bunker   Calif…            2.5     39.2     -120.
##  3 2017-01-03 00:00:00 PP Ridgetop Unit Cal F…           10       38.8     -123.
##  4 2017-01-03 00:00:00 FIA Forest hill  Priva…           10       39.0     -120.
##  5 2017-01-03 00:00:00 28A              US Fo…           51       37.7     -119.
##  6 2017-01-04 00:00:00 Plum 52          US Fo…            1       39.5     -121.
##  7 2017-01-04 00:00:00 A-01             US Fo…            5       37.8     -119.
##  8 2017-01-04 00:00:00 PP Ridgetop Unit Cal F…           10       38.8     -123.
##  9 2017-01-04 00:00:00 5                US Fo…           59       37.6     -119.
## 10 2017-01-05 00:00:00 Big Chico Creek… CSU C…            0.1     39.8     -122.
## # ℹ 10,640 more rows
## # ℹ 3 more variables: `Fuel Type` <chr>, `Total Tons` <chr>, `Burn Type` <chr>

Fix column names

xy_new <- xy_new %>% 
  rename(Burn_Date = `Burn Date`) %>% 
  rename(Burn_Unit = `Burn Unit`) %>% 
  rename(Acres_Burned = `Acres Burned`) %>% 
  rename(Fuel_Type = `Fuel Type`) %>% 
  rename(Total_Tons = `Total Tons`) %>% 
  rename(Burn_Type = `Burn Type`) 

Fix dates that should be negative

xy_new <- xy_new %>% 
  mutate(Longitude = ifelse(Longitude > 0, Longitude*(-1), Longitude))

Add year column

xy_new <- xy_new %>% 
  mutate(Year = year(Burn_Date))

Remove duplicates

xy_new <- xy_new %>% 
  distinct(Burn_Date, Burn_Unit, Acres_Burned, .keep_all = T)

Make spatial and save

2019

Filter to 2019

xy_new_2019 <- xy_new %>% 
  filter(Year == 2019)
xy_new_2019 %>% head()
## # A tibble: 6 × 10
##   Burn_Date           Burn_Unit Agency Acres_Burned Latitude Longitude Fuel_Type
##   <dttm>              <chr>     <chr>         <dbl>    <dbl>     <dbl> <chr>    
## 1 2019-01-01 00:00:00 Forest M… Leoni…          1       38.6     -121. Slash    
## 2 2019-01-01 00:00:00 Leoni Me… Leoni…          2       38.6     -121. Natural  
## 3 2019-01-02 00:00:00 Foresthi… Calif…          0.1     39.0     -121. Brush    
## 4 2019-01-02 00:00:00 Forest M… Leoni…          1       38.6     -121. Slash    
## 5 2019-01-02 00:00:00 AIRPERSON Sierr…          2       38.9     -121. Slash    
## 6 2019-01-02 00:00:00 Wyntoon … Wynto…          2       41.2     -122. UNK      
## # ℹ 3 more variables: Total_Tons <chr>, Burn_Type <chr>, Year <dbl>

Make spatial

sf_new_2019 <- st_as_sf(xy_new_2019, coords = c("Longitude", "Latitude"), crs = 4326)

Export to shapefile

st_write(sf_new_2019, "shapefiles", "PFIRS_2019_pulled2023", driver = "ESRI Shapefile", delete_layer = T)
## Deleting layer `PFIRS_2019_pulled2023' using driver `ESRI Shapefile'
## Writing layer `PFIRS_2019_pulled2023' to data source 
##   `shapefiles' using driver `ESRI Shapefile'
## Writing 1288 features with 8 fields and geometry type Point.

All years

Make spatial

sf_new <- st_as_sf(xy_new, coords = c("Longitude", "Latitude"), crs = 4326)

Export to shapefile

st_write(sf_new, "shapefiles", "PFIRS_2017-2022_pulled2023", driver = "ESRI Shapefile", delete_layer = T)
## Deleting layer `PFIRS_2017-2022_pulled2023' using driver `ESRI Shapefile'
## Writing layer `PFIRS_2017-2022_pulled2023' to data source 
##   `shapefiles' using driver `ESRI Shapefile'
## Writing 10252 features with 8 fields and geometry type Point.

Save as Rdata

save(sf_new, file = "PFIRS_2017-2022_pulled2023.Rdata")

Mask out federal lands

Crop to state of CA

CA <- st_read(dsn = paste0(ref_path, "/CA boundary/ca-state-boundary/", layer = "CA_State_TIGER2016.shp"))
## Reading layer `CA_State_TIGER2016' from data source 
##   `C:\Users\ctubbesi\OneDrive - California Air Resources Board\Documents\Reference data\CA boundary\ca-state-boundary\CA_State_TIGER2016.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13857270 ymin: 3832931 xmax: -12705030 ymax: 5162404
## Projected CRS: WGS 84 / Pseudo-Mercator
CA <- st_transform(CA, st_crs(sf_new))
sf_new <- st_intersection(sf_new, CA)
## Warning in st_is_longlat(x): bounding box has potentially an invalid value
## range for longlat data

## Warning in st_is_longlat(x): bounding box has potentially an invalid value
## range for longlat data
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
load(file = "Rdata/NF.Rdata")
NF <- st_transform(NF, st_crs(sf_new))
sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
sf_new_noNF <- st_difference(sf_new, NF)
## although coordinates are longitude/latitude, st_difference assumes that they
## are planar
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
mapview::mapview(list(NF, sf_new_noNF), col.regions=list("red","blue"),col=list("red","blue"))

Save as Rdata

save(sf_new_noNF, file = "Rdata/PFIRS_2017-2022_pulled2023_noNF.Rdata")